Discrete kink dynamics in hydrogen-bonded chains I: The one-component model 
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We study topological solitary waves (kinks and antikinks) in a nonlinear one-dimensional Klein- 
Gordon chain with the on-site potential of a double-Morse type. This chain is used to describe 
the collective proton dynamics in quasi-one-dimensional networks of hydrogen bonds, where the 
on-site potential plays role of the proton potential in the hydrogen bond. The system supports a 
rich variety of stationary kink solutions with different symmetry properties. We study the stability 
and bifurcation structure of all these stationary kink states. An exactly solvable model with a 
piecewise "parabola-constant" approximation of the double-Morse potential is suggested and studied 
analytically. The dependence of the Peierls-Nabarro potential on the system parameters is studied. 
Discrete travelling- wave solutions of a narrow permanent profile are shown to exist, depending on 
the anharmonicity of the Morse potential and the cooperativity of the hydrogen bond (the coupling 
constant of the interaction between nearest-neighbor protons). 
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I. INTRODUCTION 



Hydrogen bonds (H-bonds) play a crucial role in the 
structure and the dynamics in a whole variety of systems 
ranging from ferroelectrics to biomolecules. They are of 
central importance in biology, when reactions are con- 
sidered at molecular level. In bioenergetics, they appear 
even more crucial because they enable transfers of pro- 
tons from one molecule to another one in networks or 
chains formed via hydrogen bonding [Q . 

More specifically, H-bonds [|| are interactions link- 
ing two molecules or ions, for example, O, N, F, and 
CI atoms, or in general any pair of hydroxyl groups, 
which may be denoted by X, via a hydrogen ion (proton) 
H + , forming a hydrogen-bonded (HB) bridge X-H- • X 
as shown schematically in Fig. [|. The ion to which the 
proton in this bridge (H-bond) is more tightly linked is 
called the hydrogen donor, whereas the other ion is the 
hydrogen acceptor. More precisely, the proton is coupled 
to each X~ ion through a pair ion-proton interaction po- 
tential of the standard type (Morse, Lennard- Jones, etc.) 
with an equilibrium distance tq, which necessarily has a 
finite dissociation energy as the X- • H distance tends to 
infinity. 

As usual, the total potential for the HB proton is of a 
double-well shape, but this can occur only if the motion 
of the heavy ions along the H-bond is appropriately con- 
strained, so that either (i) a sufficiently strong interac- 
tion (repulsion) between the X~ ions that does not allow 
the ions to get closer each other a distance less or equal 
to 2ro or (ii) a periodic substrate potential with period 
exceeding 2ro is additionally involved. In this way, a one- 
dimensional network of hydrogen bonds can be formed as 
a diatomic chain of alternating heavy (ion) and light (pro- 
ton) particles coupled nonlinearly (e.g., via a Morse- type 
potential, like in Ref. whereas the second- neighbor 
(ion-ion and proton-proton) interactions are involved in 



the harmonic approximation. Under certain conditions 
on the ion-ion coupling discussed in Ref. [Q, the proton 
in each H-bond of the HB chain can be found in two equi- 
librium positions separated by a potential barrier, so that 
the two degenerate ground states of the chain • • • X-H 
• • • X-H • • • X-H • • • and • • • H-X • • • H-X • • • H-X • • • 
are assumed to exist. Another important property, more 
specific for biological systems, is that the height of the 
potential barrier crucially depends on the distance be- 
tween adjacent X~ ions. Using these properties as main 
features of HB chains, a number of one-dimensional two- 
sublattice models, whose dynamic behavior is governed 
by the soliton theory j|, has been suggested and stud- 
ied extensively. These soliton-like theories are based on 
the well-known cooperativity of the hydrogen bonding, 
simply defined through the coupling of protons in the 
nearest-neighbor hydrogen bridges of the chain. 

Since the HB chain is a diatomic lattice, the mecha- 
nism of hydrogen bonding involves two types of particle 
displacements. Let Q n and q n be the displacements of the 
heavy ion and the proton in the nth unit cell of the lat- 
tice from their equilibrium positions, at which one of the 
two ground states of the chain is realized, respectively. 
These displacements are labeled according to the se- 
quence {...,Q n - 1 ,q n -- L ,Q n ,q n ,Q n+ i,q n+1 ,...}. Then 
the general and the most simple model for the proton 
transfers in such a diatomic chain can be given through 
the two-sublattice Hamiltonian, consisting of two parts 



H = H + H it 



(1) 



The first part 



IfQn + ~2 { ^ n + 1 ~ qn > +£oV(u n ,p n ) 



(2) 



1 



with 

U n — q n - - (Q n + Qn+i) and p n = Q n+ \ — Q n , (3) 

describes the proton kinetic energy, the nearest-neighbor 
proton-proton interaction, and the intrabond proton en- 
ergy that depends on the displacements of the protons 
from the midpoints in the H-bonds u n 's and the rel- 
ative distances between the nearest- neighbor ions p„'s, 
whereas the second (pure heavy-ion) part 



\MQl + l -K lon pl + ~K sub Q 2 n 



(4) 



describes the kinetic energy of the X~ ions, the coupling 
energy between the nearest-neighbor ions, and the inter- 
action energy of the X~ ions with a possible substrate 
(e.g., formed by the walls of a pore crossing a mem- 
brane). Here the overdot denotes differentiation on time 
t. The proton and ion masses are denoted by m p and 
M, respectively. Similarly, K pi K ion , and K su i> stand 
for the stiffness constants of the interaction between the 
nearest-neighbor protons, the nearest-neighbor ions, and 
the chain ions with the substrate, respectively. It is 
important that the intrabond proton energy is given in 
terms of a general double- well potential V(u, p) as a func- 
tion of two variables: u, the proton displacement from the 
middle of the hydrogen bond, and p, the relative ion dis- 
placement. If additionally this (dimensionless) function 
is normalized according to the relations V(0, 0) = 1 and 
V(±a, 0) = 0, where ±a are the positions of minima in 
this function, then £q is the barrier height of the proton 
potential in the H-bond. When the heavy ions are dis- 
placed from their equilibria, this potential is deformed, 
with its barrier top moved together with the ions. 

There have been numerous studies j| of soliton solu- 
tions to the equations of motion governed by the Hamilto- 
nian (l)-(4) including also one-component models, where 
the heavy ions X~ are assumed to be fixed p-plj. All 
these studies refer to the continuum limit, which pre- 
sumes the existence of a sufficiently effective cooperativ- 
ity of hydrogen bonding or, in other words, the inter- 
bond proton-proton coupling K p is required to be strong 
enough. However, according to the ab initio calculations 
of the proton-proton interaction in realistic HB chains 
by Godzik @, K p ~ 41 Kcal/mol A 2 . This magnitude 
appears not to be sufficient for a free propagation of the 
ionic defects along the HB chain with realistic values of 
the potential barrier height eq. The reduction of this 
barrier on the basis of the two-component modelling was 
also shown to be not enough to provide a free soliton 
regime and the ionic defects in HB chains with realis- 
tic parameter values [ p"3| were shown to be very narrow 
objects @. 

On the other hand, a rapidly increasing number of 
publications over last years (since the pioneering work of 
Peyrard and Kruskal |15j ) have demonstrated significant 
differences in the behavior of soliton solutions treated 



in the continuum limit and their spatially discrete rela- 
tives jl(|. The discrete versions of the partial differen- 
tial equations brings about a number of critically impor- 
tant modifications to the dynamics. The moving kinks 
of the continuum theories become propagating structures 
that decelerate by emitting radiation as they traverse the 
lattice sites. This ultimately brakes the structures and 
brings them to rest, or "pins" them. 

The above results has been obtained for the conven- 
tional models such as the discrete sine-Gordon and (j) 4 
chains. For more general class of models, some inter- 
esting and intriguing results has been obtained. Thus, 
it has been shown that the shape of the on-site (in our 
case, the intrabond proton energy) potential is a factor 
of particular importance for modelling soliton motion in 
physical systems. To study this effect, Peyrard and Re- 
moissenet fl7j| have introduced a modified sine-Gordon 
system, where the shape of the on-site potential differs 
sufficiently from the sin-function. They found that if the 
barrier between the potential wells is flat enough, the 
Peirels-Nabarro (PN) barrier does not decrease mono- 
tonically with the coupling constant, as in the ordinary 
discrete sine- Gordon chain. It decreases with oscillations, 
so that the PN barrier experiences dips, where it lowers 
by the order of magnitude. Later |Q, it was found that 
if the PN barrier decreases nonmonotonically, there ex- 
ist certain velocities, at which even very discrete kinks 
propagate with constant shape and velocity. Everywhere 
in between these velocities, there exist kinks with oscilla- 
tory asymptotics (nanopterons). Approaching the prob- 
lem from another side, Schmidt in [l9| ] has constructed 
a Klein-Gordon model that allows an exact moving kink 
solution of the form tanh(n— vt) for some specific value of 
velocity v. Furthermore, Flach and coauthors have 
shown that for this model the PN barrier is nonzero 
(when « / 0). It was also shown that kinks of the 
discrete sine-Gordon equation with topological charges 
greater then one exhibit some features, similar to those, 
described above, including free propagation at some se- 
lected velocities (see Refs. |l^^,|l^]). Note that if we 
step out from the Klein- Gordon class of discrete mod- 
els (for instance, by introducing anharmonicity into the 
interparticle interaction), some new phenomena related 
to the kink mobility can appear, but this is out of the 
scope of the present paper. In this context, a few pa- 
pers (see Refs. 21 23] and others therein) should also be 



mentioned. 

Thus, owing to the importance of discreteness effects 
in the kink dynamics, it would be of big interest to ap- 
ply these findings first for the one-component model of 
proton transport in HB chains and then for the two- 
component model given by the Hamiltonian ( |l|)-(4|). In 
this context, as found by Duan and Scheiner J24]25| ], a 
pair of Morse functions, placed tail-to-tail so as to al- 
low for the approach of the proton towards the acceptor 
while it is departing from the donor (see Fig. 1), provides 
the best framework for reproducing their potentials ob- 
tained from ab initio calculations. It is important that 
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the Morse-type functions contain parameters with clear 
physical meaning, which vary little from one H-bond to 
the next one in HB systems. 

The aim of this paper is to investigate the properties of 
the one-dimensional Klein-Gordon chain with the on-site 
(intrabond) potential of the double-Morse type. We are 
going to find both stationary and dynamic (moving) kink 
solutions and to show that kinks can be mobile even if 
being very narrow. 

The paper is organized as follows. In the next sec- 
tion, we present the Hamiltonian and the equations of 
motion for the one-component model. In Sec. Ill, we 
study the properties of the stationary kink solutions. In 
Sec. IV, the Peierls-Nabarro potential for the kinks is 
investigated. Section V is devoted to the studies of kink 
mobility. Conclusions are given in Sec. VI. 



II. THE DOUBLE-MORSE PROTON POTENTIAL 

In the limit, when the heavy ions are fixed (immobile) 
at a same distance I forming a uniform lattice, we deal 
only with the first part of the Hamiltonian (g), where 
In — u n- In what follows we adopt the dimensionless 
description, where for the dimensionless proton displace- 
ment u n /l we keep the same notation u n and use the 
time unit t a — 1/ ^so/nip. In these dimensionless vari- 

reads 



l/\/foJm 
ables, the Hamiltonian (I 



;("■ 



n+l 



u n ) 2 + V(«n) 



(5) 



Here and in what follows, the overdot denotes the differ- 
entiation with respect to r = t/t and k = K p l 2 /e is the 
dimensionless proton-proton coupling constant. 

As described in Introduction and illustrated by Fig. 1, 
the intrabond proton potential V(u) can be formed as a 
result of superposition of two pair ion-proton interaction 
potentials placed tail-to-tail. According to the ab initio 
studies of Duan and Scheiner 24 2j| , these potentials are 
preferred to be chosen of the Morse type. As a result, the 
potential V(u) is a symmetric double- well function with 
minima at u = ±a and a maximum at u = y,EGLp7| : 



V(u) 



a — cosh (f3u) 
a-1 



a = cosh(/3a) 



(6) 



The inequality a > 1 ensures the double-well form of 
the function (^). Throughout the paper we take a = 
0.25. The potential (||) is normalized so that the barrier 
height always equals unity. Its shape for different values 
of f3 is shown in Fig. |J As can be seen from this figure, 
the parameter (3 determines the curvature/flatness of the 
barrier and its shape strongly depends on this parameter. 
We assume its values to range over the whole half-axis 
< (3 < oo. For small f3, the barrier is rather narrow 
being the limiting case of the <f> model, i.e., 



lim V(u) = (l-u 2 /a 2 ) 2 . 



(7) 



Increase of f3 makes the barrier more flat and the wells 
more narrow, so that the other limit is 



lim V(u) 



-oo < u < —a, 

u = ±a, 
—a < u < a, 
a < u < oo. 



(8) 



The corresponding equation of motion is the well 
known discrete nonlinear Klein-Gordon equation: 



in = k (u n+ i - 2u n + u n -i) - V'(u n ), n = 0,±l, . 



(9) 



Here and in what follows the prime denotes differentia- 
tion of a function with respect to its argument. 

Before embarking on a complete analysis of the discrete 
equation ([}]), we calculate the dispersion law of small- 
amplitude waves around one of the two ground states. 
As a result, this law is given by the following equation: 



+ 2k(1 — cos q), 



UJ 



1 



a — 1 



V2[3 



tanh (Pa/2) 



(10) 



The gap of the spectrum depends on the parameter f3; 
for large values of [3, it increases linearly with (3. 



III. KINK STATES, THEIR STABILITY AND 
BIFURCATIONS 

In this section, we start from the classification of possi- 
ble stationary (anti)kink states. To compute these states, 
we have used both the conjugate gradients method for 
minimization of the stationary part (u n = 0) in the 
Hamiltonian (^|) and the Newton iteration method for 
solving the time independent nonlinear set of equations 
that originates from Eqs. (ph: 



n+l 



2u n + Un-i) = V'(u n ), n = 0,±1, 



(11) 



These equations can be rewritten as a two-dimensional 
map: 



Pn+x = u n , 
-i 



u n +i = K V'(u n ) + 2u n - p n , n = 0, 1, . 



(12) 



In general, maps of this type are chaotic. However, an on- 
site potential, for which the map ( |l2| ) is integrable, has 
been found in Ref. p8[ . The previous knowledge on the 
kink solutions in the most popular cases of the discrete 
nonlinear Klein-Gordon lattices such as the sine-Gordon 
or 4 chains implies the existence of only two stationary 
kink states. These states possess the inversion symmetry 
with respect to the center of the kink, being monotonic 
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functions on the lattice. They connect two hyperbolic 
fixed points of the map (—a, —a) and (a, a), which are 
the ground states of the chain. 

Adapted to our case of a HB chain with the num- 
bering of ions and protons according to the sequence 

.}, one of these 

stationary states, which has its center positioned at a 
heavy ion (call it an ion-centered kink/antikink), say with 
a number no, is dynamically stable, whereas the other, 
with its center positioned in the middle of a noth H-bond, 
i.e., in between the n th and the (n + l)th ions (call it a 
bond- or proton-centered kink/antikink), is dynamically 
unstable. The symmetry of the ion-centered (on the noth 
ion, in between H-bonds no — 1 and no) kink/antikink is 
defined by the relations 



*n +n— 1) 



n = 0,±1, 



(13) 



whereas for the proton-centered (in the middle of 
the noth H-bond, in between ions no and no + 1) 
kink/antikink, the symmetry relation is given by 



l n +n ) 



n = 0,±1, 



(14) 



The solutions of these types certainly exist also in our 
model, as illustrated by Fig. [| However, their stability 
properties appear to be much more complicated and they 
are discussed below. 

In general, a deformation of the barrier shape in the po- 
tential V(u) leads to a more rich family of kink/antikink 
solutions. This has been observed in the previous studies 
0,^8|. The first feature caused by the deformation of 
the proton potential (|^) with increase of (3 is the phe- 
nomenon of stability switching, according to which the 
two types of kink solutions with inversion symmetry de- 
fined by Eqs. ( |l3| ) and (|lj) switch their stability, while 
varying the system parameters. The second one is the 
appearance of new types of kink solutions. For the first 
of these types the symmetry relations ([l3] ) and (|lj) are 
not valid anymore, whereas the other one is symmetric 
with a zigzag-like profile at their center, but still has 
monotonic asymptotics as |n| — > oo. 



A. An exactly solvable limit 

To understand better the effect of stability switching, 
it is instructive to consider the limit (3 — > oo, resulting 
in a similar potential behavior as studied in Ref . |l8| . In 
this limit, the double-Morse potential takes the form [see 
Eq. (||)] , for which the system of equations of motion ( p"l| ) 
becomes exactly solvable because the particles (protons) 
can appear only either in the wells or on the flat region of 
the proton potential V(u). There exists an infinite, but 
a countable set of the stationary kink solutions with an 
arbitrary number of particles, m = 0, 1, . . ., lying on the 
barrier. This number uniquely defines a kink/antikink 
solution. The set of these solutions can be written as 




-oo < n < no 



n 



2 

n 



- 1 < n < n 4 
f < n < oo, 



(15) 



for the kink centered on the noth ion (rn = 0, 2, . . .) and 



m+l ' 

a, 



"o 



-oo < n < no 

m+l 
2 

n a + < n < oo, 



m+l 

2 : , 
m+l 

2 



m+l ^ ^ i m+l / -i r\ 

< n < n -I ^— , (16) 



for the kink centered on the n th H-bond (m =1,3,.. .). 
Similarly, the analytical expressions for the antikink so- 
lutions can be obtained. These stationary solutions can 
be found, using either energy arguments based on the 
Hamiltonian (^|) or directly from the equations of motion 
(liTI), where 







lim V'{u) = 



—oo, 
0, 
oo, 



— oo < u < —a, 
—a < u < a, 
a < u < oo. 



(17) 



The energy of both these (anti)kink solutions in the 
limit (3 — + oo is easily calculated and for any integer 
rn = 0, 1, . . . , it reads 



E m = E m (K) = m + 2na 2 /(rn +1). 



(18) 



As illustrated by Fig. ^, where the linear dependences of 
the energy E m on k are plotted for different m's that the 
crossings of these dependencies occur at some values of 
the coupling parameter n. 

The most interesting points in Fig. ^ are the crossings 
for the states with the energies that correspond to ad- 
jacent m's, i.e., m = and 1, m = 1 and 2, and so on 
because they occur at lowest energies. Thus, the cross- 
ings when the kink state with m particles on the barrier 
is transformed into the state with m + l particles on it, 
occur at the following values of k: 



M 

v m+l 



= (m+ l)(m + 2)/2a 2 , m = 0, 1, . 



(19) 



Therefore, depending on the strength of the proton- 
proton coupling k, the proton-centered or the ion- 
centered kink can reach a global minimum of the energy 
E m > (except for the ground states, when E = 0). In 
addition, it is interesting to notice that at the values of 
k given by Eqs. (|l9|), the inter-bond and the intrabond 
energies are equal exactly each other. 



B. Numerical results for finite /3's 

Now let us investigate how the properties of the kink 
solutions found in the exactly solvable limit (3 — > oo 
change, when (3 take finite values. Thus, changing (3 
allows us to explore the whole set of scenarios, starting 
from the </> 4 limit and finishing with the limit (3 — * oo. 
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We start to compute the kink solutions from the anticon- 
tinuous limit (k = 0), taking the solutions of the exactly 
solvable limit as an initial guess for the Newton iteration 
method. Then, we increase k and check how the kink 
profiles behave. In Fig. |^, the energy dependence on the 
coupling parameter k is plotted for (3 = 10. First, let 
us focus on the behavior of the solutions with the lowest 
energies (namely those which correspond to to = and 
m = 1 in the exactly solvable limit). 

In Fig. |E|, curve 1 corresponds to the ion-centered kink 
with m = 0, the symmetry of which is given by Eqs. (13), 
whereas curve 2 to the proton-centered kink with m = 1, 
the symmetry of which is given by Eqs. (|l4|). Contrary 
to the previous knowledge on the stability properties of 
the kink solutions for the discrete sine-Gordon and </> 4 
models, where the (anti)kinks of the symmetry (^) are 
always dynamically stable, while the (anti)kinks of the 
symmetry ( |l4| ) are always dynamically unstable, an in- 
terchange of stability is observed for the proton potential 
(^J) with finite /3's as K varies. Thus, one can see that at 
a certain value of the coupling parameter, k ~ 24.5, the 
energies of the both types of symmetric kinks coincide 
and after passing this critical point, the proton-centered 
kink appears to be stable, while the ion-centered kink is 
unstable. For higher k's, several more interchanges of 
stability take place, with the energy difference between 
the sequential kink states that decreases with the growth 
of k. These transitions of stability take place smoothly 
all the way up to the continuum limit. We refer to these 
transitions as to stability switchings. Thus, the solutions 
with to > 1, which were clearly separated from the m = 
and the m = 1 solutions in the exactly solvable limit 
(3 — ► oo, appear to be smoothly connected with them. 
In other words, while the coupling k increases, the par- 
ticles slowly "climb" on the barrier, so that there is no 
abrupt transition from the state with to = to the state 
with to = 2 or, further, to the states with m = 4, 6, . . .. 
The same can be concluded about the kinks with to be- 
ing odd. In this respect, the system still shows similarity 
with the (j> model. 

If we focus more carefully on the behavior of the system 
in the vicinity of the points, where the energies of proton- 
centered and ion-centered kinks become equal, we find 
that a new type of kinks appears. These kinks shown 
in Fig. U do not exhibit any of the symmetries defined 
by Eqs. ([13]) and ([li]), but they are doubly degenerate 
related to each other by the inversion with respect to 
the crossing point between the line u n = and the line 
connecting two central particles of the kink [n — 48 and 
n = 49 in Fig. g(a)]. The energy of these kink states 
with broken symmetry is always larger than the energy of 
the symmetric (ion-centered and proton-centered) kinks. 
The asymmetric (anti)kinks are always linearly unstable. 
In Fig. ||, their energy is shown by curve 3. 

The kinks of the last type are shown in Fig. |5| by curves 
4 and 5. This is what happens to the solutions with to > 
1 obtained in the limit /3 — * oo, when we are moving from 
the anticontinuous limit (n = 0). Instead of attaining a 



regular monotonic form, the kink profiles with several 
particles on the barrier develop a zigzag-like structure in 
their centers, as demonstrated by Fig. 0. For finite /3's the 
zigzag-like solutions are linearly unstable and therefore 
we do not study them here in more detail. 

Instead, we focus on the behavior of the monotonic 
kink solutions with increase of the coupling k in the vicin- 
ity of the stability switchings. To understand this effect 
better, in Fig. ||, we have plotted the position of the 
(iV/2)th proton (N is the total number of H-bonds in 
the chain) as a function of the coupling n. Here the se- 
quence of pitchfork bifurcations is clearly seen. Curve 
1 corresponds to the kink centered in the middle of the 
(7V/2)th bond, i.e., when the (A/2)th proton is always a 
central particle of the kink with ujv/2 = 0. The displace- 
ments of the (iV/2)th proton for the ion-centered kinks 
positioned on the (N/2 + l)th and the (N/2)th ions are 
shown by curves 2 and 3, respectively. When increasing 
k, the (iV/2)th proton moves slowly out of the well. At 
a certain value of k, more specifically, at k ~ 22.5, the 
pitchfork bifurcation of the proton-centered kink takes 
place. This configuration retains its stability and two 
new solutions (both linearly unstable) appear. These are 
precisely those asymmetric kinks (see curves 4 and 5 in 
Fig. ||), with their shape shown in Fig. [| At the be- 
ginning, they look like slightly distorted proton-centered 
kinks, but with the growth of k, they change more and 
more towards the ion-centered configuration. Eventually, 
the second pitchfork bifurcation takes place at k ~ 26.2. 
The asymmetric kinks join the ion-centered kinks (junc- 
tion of curves 3 and 4, and curves 2 and 5) and the ion- 
centered configuration loses its stability. Now one can 
clearly see that two identical kinks, shifted by one lat- 
tice spacing with respect to each other, are connected 
via such a bifurcation sequence. Thus, this pitchfork bi- 
furcation is nothing but a transition of the kink from one 
position to another position, one lattice period forward 
or backwards. This cascade of bifurcations can be contin- 
ued further up or down in u n 's or, in other words, two or 
more sites backwards or forward. Similar bifurcation sce- 
nario for discrete breathers in the ac-driven and damped 
Klein-Gordon lattice has been reported in p9| . 

For higher values of /3, the effect of stability switch- 
ings exists, being more pronounced because the switch- 
ings start at smaller n and take place more frequently 
(see Fig. ^). Another feature that appears from the ex- 
actly solvable limit (3 — > oo is as follows. Curve 1 of 
Fig. H corresponds to the ion-centered kink and curve 2 
to the proton-centered kink. They cross each other at 
k ~ 17.4, where the ion-centered kink loses its stability, 
disappearing shortly at k ~ 28.6. Here the same pitch- 
fork bifurcation scenario takes place, but the distance be- 
tween the first and the second bifurcations is much larger 
than in the [3 = 10 case. In the meanwhile, a bit earlier, 
at k ~ 23.4, a new family of ion-centered kinks appears 
(curve 4). This curve corresponds to the ion-centered 
kinks with two protons on the barrier. Thus, one can 
observe the coexistence of two different kink solutions 
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with the same type of symmetry (for more details see 
the upper inset of Fig. . This coexistence takes place 
on a rather narrow interval of k and one of the coexist- 
ing kinks is unstable, but still this phenomenon clearly 
originates from the limit (3 — ► oo. In Fig. we show 
the shape of two coexisting kinks. One of them (more 
narrow) corresponds to curve 1 of Fig. ^ and the second 
one, which is more broad, corresponds to curve 4 of this 
figure. These kinks have zero and two protons on the bar- 
rier, respectively. When k increases further, the stability 
switchings occur between curves 2 and 4 under the same 
scenario as before for (3—10 (see for details the lower 
inset of Fig. ^|). Zigzag-like kinks are also presented in 
this shown by curves from 6 to 8. 

When increasing (3 even more, the coexistence of dif- 
ferent kinks with the same symmetry is even more pro- 
nounced. We have checked the case of (3 = 50 and discov- 
ered for this value that several cases of this coexistence 
for different kinks take place and they are much more pro- 
nounced. Thus, this case is rather close to the exactly 
solvable limit j3 — *■ oo. 



C. Analytical approximation for finding kink 
solutions 

One can use a simplification of the potential (|^) in or- 
der to obtain analytically an exact (anti)kink solution. 
To this end, we approximate both the wells of the poten- 
tial by parabolas connected by the constant equal to the 
barrier height as follows: 



V(u) 



H/2){u- 
1, 



a) 2 , 
a) 2 , 



— oo < u < —6, 
-b < u < b, 
b < u < oo, 



(20) 



where b = a — \f2jujo- Schematic description of 
this "parabola-constant" approximation is presented in 
Fig. |ll| by thick solid lines. The thin solid lines show the 
original potential (^) with (3 — 20. The approximation is 
expected to work well, when the barrier is flat enough, 
i.e., when loq 3> \[2/a. Within this approximation, we 
are able to solve the problem of finding stationary kink 
solutions analytically. 

Let to — 0, 1, . . . , be the number of protons on the bar- 
rier of the potential ( |20| ) . Then the discrete kink profiles 
are given by 



+ y4 TO g A (" _,l O+ m / 2 + 1 ) 



if 



= < 



oo < n < no — m/2 — 1, 
(n-n + l/2)D m 
if no — m/2 — 1 < n < no + m/2, 
a — A m e~ x ( n ~ n °~ m '/ 2 ) 
if no + m/2 < n < oo, 



(21) 



for the kink centered on the ngth ion (m = 0, 2, . . .) and 



if n 



if — oo < n < no — (m + l)/2, 

(n - no)D m , , 

(m + l)/2<n<n Q + (m + l)/2, 1 ' 
a — A e~^ n ~ n o-(m+i)/2] 
if no + (m + l)/2 < n < oo, 



for the kink solution centered on the noth H-bond (m = 
1,3,.. .). Here A is a "localization" parameter that mea- 
sures the transition width in the (anti)kink profile be- 
tween the uniform distribution of protons on the barrier 
and the (anti)kink asymptotics u n — * ±o. It is given by 
a positive root of the equation 



cosh A = 1 + Uq/2k. 



(23) 



The other two parameters, the amplitude A and the uni- 
form distance between the nearest-neighbor protons on 
the flat D can be expressed through the localization pa- 
rameter A as 
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D m = 



2a 



(rn + l)e A — to + 1 ' 
2ae x 

(to + l)e A — to + 1 



(24) 



As follows from Eqs. (|T^) and (^3|), the limit A — > oo 
(when Wq/k — > oo) is more general than the limit (3 — > oo 
because A contains both (3 and k. Therefore one can 
check that Eqs. (|2l])-(24) are reduced to the stationary 
kink solution given by Eqs. (|l^) and ( |l6| ) as A — > oo. In 
particular, the amplitude A m and distance between the 
protons on the barrier D m tend to zero and 2a/ (m + 1), 
respectively. 

Using Eqs. ([H) 
ergy of both the kin. 



one can easily compute the en- 
c configurations: 



E,, 



m + 2na 



2 tanh(A/2) 
1 + TOtanh(A/2) ' 



(25) 



Similarly, this expression is also transformed to the en- 
ergy ( |l8| ) as A — > oo. 

Now we investigate the behavior of the energy differ- 
ence 



AE m+1 (n, (3) = E m+1 - E m , m = 0, 1, . 



(26) 



We find that this difference as a function of k has a num- 
ber of zeroes and these zeroes depend on the parameter 
(3. Thus, within our approximation one can predict the 
effect of switching of the stable and unstable kink config- 
urations. In Table |, we show the values of the coupling 
parameter k, for which the first [i.e., when k — K\ and 
to = 0; see also Eq. ([H)])] switching of the kink stability 
takes place. 

TABLE I. Comparison of numerically and analytically cal- 
culated values of k for which the first (m = 0) stability switch- 
ing occurs. 



G 






Ki, numerical 


Ki, analytical 


5 


64.10 


18.0 


10 


24.547 


17.069 


20 


17.383 


16.329 


50 


16.172 


16.052 


oo 


16.0 


16.0 



We see that the approximation works fairly well when 
(3 is rather large and the barrier between the wells is close 
to being completely flat. It is improving with increase of 
(3 and in the limit (3 — > oo, our approximation coincides 
with the exact result, shown above. Thus, switchings of 
the stability of kink states with different symmetries is a 
generic effect that does not depend on a specific model, 
but on the properties of the on-site potential. Another 
piecewise approximation of the proton potential V(u) for 
a HB chain has been constructed earlier by Weiner and 
Askar p5| , using alternatively inversed parabolas. Our 
potential approximation ( ppf ) seems to be more appro- 
priate for the studies in detail of flatness effects of the 
on-site (intrabond) potential and more close to the real- 
istic double-Morse potential (|J) if (3 is not so small. 

Note that the zigzag-like kink profiles obtained above 
numerically shown, e.g., in Fig. 7 can also be given an- 
alytically within the approximation (p0|). Indeed, the 
ion-centered kink shown in Fig. 7(a) is described by 

-a + B e x( * n+2 \ n= -2,-3,..., 

u = = £ Q , (27) 

a- Boe-^"-^, n= 1,2,..., 

whereas the proton-centered profile illustrated by 
Fig. 7(b) is given by 

-a + B x e x{ - n+2 \ n = -2,-3,..., 

u = 0, mi = -u-i = £i, (28) 
a- Sie- A (™- 1 ), ti = 1,2,..., 

where A is given by Eq. ( p3|) and 

" X e~ A - m/2 - 1 



B m — 2a- 



(e x — m + 1) 

3 — p x — 2p~ x 

tm=a— — , m = 0,l. (29) 

e A — m + 1 

Using the last equations, the kink energy can be calcu- 
lated in a similar way as the energy (|l8|). In the limit 
A — > oo, the energy of the zigzag-like kinks, profiles of 
which are illustrated by Fig. 7, become 

E Q = Qua 2 and E x = 1 + 5na 2 . (30) 

These expressions clearly show that for this type of kinks 
their energies also can coincide for certain values of cou- 
pling k if f3 (and, consequently, A) is large enough. 



D. Elementary excitations on the kink background 

Let Un^ be a stationary kink solution of Eq. (p"i"|), i.e., 
a fixed point of the map (fl2"|). We are interested in 
the properties of small-amplitude excitations on the kink 
background. Linearizing Eq. (||) around the stationary 
kink solution according to 

u n = u^+A n e m \ (31) 
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we arrive at the eigenvalue problem 

LA = AA, A = {. . . , A n -i, A n , A n+ i, . . . ,}. (32) 

Here the operator L [the Hessian of the Hamiltonian 
Eq. (|J)] acts on a vector A as 

(LA) n = -K(A n+1 - 2A n + An-x) + V n A n , (33) 

where V n — V"[un ]■ The operator L is a symmetric (so 
all eigenvalues are real) tridiagonal matrix and the spec- 
tral parameter is A = fl 2 . This eigenvalue problem can 
be treated as a quantum-mechanical problem of a par- 
ticle, trapped in a single-well spatially discrete potential 
formed by the kink. Its depth depends on the curva- 
ture/flatness of the proton potential in the middle of the 
H-bond, V"(0), and tends to uj 2 as n — > ±00. 

The eigenvalues of the problem also give information 
about the linear stability of the kink solution. If there 
exists at least one eigenvalue A = il 2 < 0, the linear ex- 
citation on the kink grows exponentially in time and the 
corresponding kink solution is linearly unstable. Oth- 
erwise, it is linearly stable. The stability of stationary 
kink solutions is determined by the system parameters 
(in our case, by the curvature/flatness parameter j3 and 
the coupling n). 

In Fig. [L2|, we depict the dependence of the eigenfre- 
quencies f2 n 's on the coupling parameter k for different 
curvatures (3. The spectrum consists of two parts. The 
first one describes the eigenfrequencies that lie above the 
potential V n and therefore determine delocalized eigen- 
vectors corresponding to the linear (phonon) spectrum 
of the lattice. The second part consists of eigenfrequen- 
cies that lie below uj 2 corresponding to spatially localized 
eigenvectors. The lowest eigenfrequency is the Goldstone 
mode, which is universally presented in all the Klein- 
Gordon models. 

The second mode corresponds to small-amplitude os- 
cillations of the kink core, being often referred to as the 
Rice mode. In the continuum cj) model, there are only 
two localized modes: the Goldstone mode and the Rice 
mode. The properties of the internal modes can be signif- 
icantly altered due to different factors such as change of 
the nature of the interparticle interaction |3l| or change 
of the shape of the on-site potential |32| . If /5 is not so 
large, the behavior of the eigenfrequencies and eigenvec- 
tors (see Fig. 13) of our system is reminiscent to that of 
the 4 model, as shown in panels (a) and (b) of Fig. |l2] . 
However, several differences occur. Thus, the Goldstone 
mode collides with the zero axis and becomes unstable 
[see Fig. Il2](a)] for the ion-centered kink. Meanwhile, for 



the proton-centered kink [see Fig. 12(b)], the Goldstone 
mode was initially unstable and became stable later on. 
The stability switchings are caused by the pitchfork bi- 
furcations as described at the beginning of this section. 
Another difference is the appearance of the new eigenfre- 
quency for the proton-centered kink with the eigenvector 
which has two nodes [see Fig. |l^(e)] . 



Thus, one can see that curvatures play an important 
role in the properties of the eigenfrequencies of our sys- 
tem. Increasing (3 further, we observe more pitchfork 
bifurcations of the Goldstone mode and the appearance 
of more localized modes. In panels (c) and (d), we show 
how the eigenfrequencies behave for /3 = 10. We observe 
more localized internal modes, some of them surviving 
in the continuum limit and some of them disappearing 
there. The shape of the corresponding eigenvectors (see 
Fig. |lj) behaves accordingly to the wavefunction shape of 
the bound states of the quantum- mechanical Schrodinger 
equation. 



IV. PROPERTIES OF THE PEIERLS-NABARRO 
BARRIER 

In general, the propagation of topological solitons 
(kinks and antikinks) in lattices are subject to their dis- 
creteness. The discreteness effects can be described by 
a spatially periodic potential with the period coinciding 
with the lattice spacing I, known as a Peierls-Nabarro 
(PN) potential. The PN potential -EpAr(n c ) is a function 
of the kink center position 



Un+l — Un-1 



1 



2(u c 



) -^Z^ n ( u n+l-Un-l), 



and the height of the PN barrier equals 



(34) 



(35) 



This difference measures the activation energy for the 
kink propagation though one lattice spacing. Here n cmax 



and 



are the positions of the kink, where the PN 



potential has its maxima or minima, respectively. The 
minimum position n c ^ m i n coincides with the kink posi- 
tion, which is a global minimum of the kink energy. Thus, 
normally n CjJTl i„ is an integer or a half-integer. 

It is a well established fact that for the </> 4 and sine- 
Gordon models, the height of the PN barrier coincides 
with the difference between the energies of the proton- 
centered and ion-centered kinks. As demonstrated above, 
the deformation of the proton (intrabond) potential V(u) 
leads to the switching of stability of these two states. 
This does not mean, however, that the PN barrier disap- 
pears, when the energy of these states coincides because 
the site-centered and proton-centered states do not rep- 
resent the states with the highest and the lowest energy 
of the PN potential 0,|l|,||. 

In Fig. [L5[ some examples of variation of the energy 
E(n c ) for different values of k are presented. Panel (a) 
corresponds to the situation, when the proton-centered 
kink is an energy minimum and the ion-centered one 
reaches a maximum. The second panel (b) demonstrates 
the case, when the first pitchfork bifurcation occurs and 
the asymmetric kink appears. Panel (c) corresponds to 
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the case, when the stability switching happens and both 
the proton- and ion-centered kink states get minima, and 
the asymmetric kinks are maxima. The period of the PN 
potential is decreased by one half. After that, as shown 
in panel (d), the proton-centered kink becomes unstable 
and the ion-centered one stable. The last panel (e) is ob- 
tained after the second pitchfork bifurcation has occured 
and the asymmetric kinks disappear. 

Dependence of the PN barrier is non-monotonic not 
only in n but also on (3. In Fig. [l6|, we show the depen- 
dence of the PN barrier on the parameter (3 for a fixed 
value of k. We observe a minimum at (3 ~ 8.3, which 
is obviously caused by the stability switching. Indeed, 
as shown by the upper inset, the minimum of the PN 
potential moves from a half-integer value to an integer 
one. Thus, the stable kink configuration changes from 
the ion-centered to the proton-centered kink. Note, that 
the height of the PN barrier does not attain zero (see the 
lower inset), but decreases by two order of magnitude. 



V. DISCRETE KINK MOBILITY 

Another interesting consequence of the deformation of 
the on-site potential barrier is the possibility of the ex- 
istence of non-oscillating travelling kinks. In the con- 
tinuum Klein-Gordon models which admit moving topo- 
logical solitons (kinks and antikinks), the domain of ad- 
missible soliton velocities s is the interval < s < Co, 
where c = \J~k, is the characteristic velocity. Thus, kinks 
in these continuum models are a one-parametric family 
of solutions with the kink velocity s as a parameter. In 
a general case of the discrete Klein-Gordon model, this 
family is reduced to a discrete set of travelling kink solu- 
tions with some velocities so, s\, . .., Sk, so = 0, Sk < cq. 
In the 4> A and the sine-Gordon models, there exists only 
sq. In general, everywhere in between s ra 's, there ex- 
ist moving kinks with oscillating asymptotics known as 
nanopterons. The existence of velocities s n 7^ has been 
shown both numerically |ll| and analytically |l9|,^0| for 
several models. This effect is due to the properties of the 
on-site (intrabond) potential V{u) and, more precisely, 
the flatness of its barrier which, in its turn, causes the 
stability switchings described in the previous sections. 

For finding non-oscillating kink solutions, we have used 
a pseudospectral method 1 33 1 . The method allows us to 



find the travelling-wave kink solutions of the type 
u n (t) = u(n — st) = u(z), 



(36) 



solving the differential equation with advanced and delay 
terms: 

sV'(z) = k[u(z + 1) - 2u{z) + u(z - 1)] - V'[u{z)]. 

(37) 

Dependence of the kink velocity on the coupling param- 
eter k is shown in Fig. 17. As follows from this figure, 



it starts (see curve 1) at the value of n, which is close 
to that when the pitchfork bifurcation takes place (see 
Fig. ||). Then the velocity grows with n. The second 
velocity dependence starts at n being close to the point 
of the second pitchfork bifurcation. 

Dependence of the kink velocity on the parameter (3 
has similar behavior, as demonstrated by Fig. |fj|. It is 
interesting to note that the value of the velocity of the 
moving kink does not decrease down to zero as n or (3 
decreases, but stops at some finite s instead. This means 
that the kinks should have sufficient kinetic energy to 
overcome the pinning effects. 

In Fig. [1^, we plot the profiles u(z) of the moving kinks 
that correspond to both two curves of Fig. [l7[ The profile 
in panel (a) corresponds to the moving kink from curve 

1, whereas the profile in panel (b) to the kink from curve 

2. In both the cases, the coupling constant was the same: 
k = 100. The second kink appears to be wider and this 
can be explained by the fact that its velocity is more 
than two times smaller. The insets of the figures show 
the velocity profiles. 

In Fig. ^0|, we investigate the behavior of moving kinks 
for different /3's. Increase of the barrier flatness is re- 
flected in change of the shape at the kink center. We 
have considered the kinks corresponding to curve 2 of 
Fig. Eh, when the coupling is fixed: k — 120. Deforma- 
tion of the kink profile can easily be seen on the velocity 
profile shown in the inset. The kink profile experiences 
deformation of its slope part, which is seen as a dip in 
the velocity profile. The dip grows with increase of (3. 
More flat the barrier becomes, more possibilities of the 
kink's profile deformation. 

These results demonstrate that the existence of a fi- 
nite set of velocities, at which the kink can move with 
constant shape is a generic effect. Its appearance is due 
to the shape of the on-site (intrabond proton) potential. 
If the proton potential has a barrier which is flat enough 
to allow the symmetry switchings (accompanied by the 
pitchfork bifurcations) , the Peierls-Nabarro barrier expe- 
riences lowering for specific values of the system param- 
eters. Thus, it becomes possible that the kinetic energy 
of the kink is sufficient to overcome the pinning forces of 
the lattice. 



VI. CONCLUSIONS 

We have studied the dynamics of the one-dimensional 
Klein-Gordon lattice with the on-site potential of the 
double-Morse type. This is a physically motivated model, 
which is the simplest one for the proton transport in a 
hydrogen-bonded chain, where the on-site potential plays 
the role of the potential for proton transfers in the hy- 
drogen bond. Therefore throughout this paper we call 
it an it intrabond proton potential. A Morse-type func- 
tion was found to offer the best combination of accuracy 
in reproducing quantum-mechanically computed poten- 
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tials [|24|,^5|. The model has two parameters, the proton- 
proton coupling k and the anharmonicity of the Morse 
potential, the curvature parameter (3. The anharmonic- 
ity parameter is responsible for the shape of the intra- 
bond potential V(u), especially on the convexity of its 
barrier. For larger (3, the barrier becomes more fiat and 
the wells become more narrow. Changing this param- 
eter, one can explore the variety of possible intrabond 
potentials, starting from the 4> 4 model (as (3 — > 0) and 
finishing with the exactly solvable limit (3 — > oo. 

Deformation of the barrier of the intrabond potential 
becomes crucial for the properties of the stationary kink 
solutions. While the </> 4 limit allows us only two types of 
stationary kinks: ion-centered and proton-centered with 
their stability properties being constant for any coupling 
k, the opposite limit shows the existence of an infinite 
countable set of stationary kink solutions. For some val- 
ues of the coupling n forming an infinite countable set, 
the states with different symmetries can have the same 
energy. In between these two limiting cases, some of the 
kink properties survive from the exactly solvable limit 
(3 — > oo. First of all, the stability properties of the kink 
solutions depend drastically on the coupling parameter 
k. With increase of k, the initially stable ion-centered 
kink becomes unstable, while the proton-centered kink 
gets stable. With further increase of k, these stability 
switchings which are, in fact, pitchfork bifurcations, go 
on several times. Another result can be seen for rather 
high /3, at least, for (3 > 20, in our calculations. It is the 
coexistence of several kink solutions of the same symme- 
try for the same k. This is a left-over from the exactly 
solvable limit and it disappears with lowering (3. An an- 
alytical approximation has been constructed to show the 
effect of the symmetry switchings analytically and to con- 
firm that the effect is not confined to a specific model, 
but has a universal nature. 

The stability switchings contribute to the non- 
monotonic behavior of the Peierls-Nabarro (PN) barrier. 
The barrier decays with increase of n, however, expe- 
riences local minima at the stability switching points, 
where its value decreases by one order of magnitude. The 
same happens for the dependence of the PN barrier on (3. 
This phenomenon, in its turn, assists the kink mobility 
and leads to the appearance of a finite set of velocities for 
which the propagation of very narrow kinks is possible. 
This is a generic effect, attributed to the shape of the on- 
site (intrabond) potential, and also is not confined to a 
specific model. Note that the PN barrier is not required 
to vanish completely (see Ref. |^(|). Simply, the barrier 
is low enough for the kinetic energy of the kink to carry 
it over the barrier. 

Thus, we have concluded that highly mobile kinks 
are possible in our model of the proton transport in 
hydrogen-bonded chains even for those proton-proton 
nearest-neighbor interactions, when the proton kinks are 
very narrow. This gives us a reason to believe that the 
soliton mechanism of proton transfers can work for phys- 
ically reasonable values of the proton-proton coupling k. 
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FIG. 1. Schematics of interactions in the hydrogen bond. 

FIG. 2. The shape of the intrabond potential V(u) given 
by Eq. (|) with a = 0.25 for (3 = 5 (curve 1), = 20 (curve 
2), and (3 = 50 (curve 3). 

FIG. 3. Profiles of monotonic symmetric kinks with /3 — 5 
and re = 30: (a) ion-centered kink and (b) proton-centered 
kink. 

FIG. 4. Dependence of the kink energy E m , m — 0, 1, 2, 3, 
and 4, given by Eq. (|l8|), in the exactly solvable limit (3 — > oo 
on the coupling parameter re. 

FIG. 5. Dependence of energy for /3=10 on coupling re 
for symmetric ion-centered kink (curve 1), symmetric pro- 
ton-centered kink (curve 2), and zigzag-like kinks (curves 4 
and 5). The inset shows more detailed behavior in the vicin- 
ity of stability switchings and curve 3 corresponds to the 
kink with asymmetric profile. Solid lines correspond to stable 
states and dashed lines to unstable ones. 

FIG. 6. Profiles of monotonic asymmetric kinks with 
j3 = 10 and k — 24. Inversion of these profiles is clearly 
seen from comparison of panels (a) and (b). 

FIG. 7. Zigzag-like kink profiles for (3 = 10 and re — 8: (a) 
ion-centered and (b) proton-centered kinks. 

FIG. 8. Dependence of the displacement of the central par- 
ticle from equilibrium (at u N / 2 = 0) on the coupling re for 
proton-centered (curve 1), ion-centered (curves 2 and 3), and 
asymmetric (curve 3) kinks with (3=10. Solid lines show sta- 
ble kinks and dashed lines unstable ones. 

FIG. 9. Dependence of the kink energy for [3=20 on the 
coupling re (see text for details). The solid line shows stable 
states and the dashed lines unstable ones. 

FIG. 10. Kink solutions corresponding to curves 1 (o) and 
4 (+) of Fig. | for re = 25. 

FIG. 11. Schematic representation of the approximate po- 
tential @. 
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FIG. 12. Dependence of the system eigenfrequencies Q n on 
the coupling parameter n for (5 — 5: (a) ion-centered and (b) 
proton-centered kinks and for f3 = 10: (c) ion-centered and 

(d) proton-centered kinks. Curves depicted by small dots cor- 
respond to cases, when Q n 's are purely imaginary and there- 
fore Im f2 n 's are plotted instead (see text for details). 

FIG. 13. Example of eigenvectors for the case with (3 — 10 
and k — 10. For the ion-centered kink: (a) the eigenvector 
of the lowest localized mode and (b) the eigenvector of the 
first excited localized mode. For the proton-centered kink: 
the eigenvectors of (c) the lowest, (d) the first, and (e) the 
second localized modes. 

FIG. 14. Example of eigenvectors for (3 = 10 and n = 30: 
(a), (b), (c), and (d) - eigenvectors of the localized modes for 
the ion-centered kink from the lowest to the highest mode; 

(e) , (f), and (g) - the same for the proton-centered kink. 

FIG. 15. Peierls-Nabarro potential for /3 = 10 and (a) 
k = 70.5, (b) k = 71.22, (c) k = 71.288, (d) k = 71.35, and 
(c) k = 72. 

FIG. 16. Dependence of the height of the PN barrier on 
the parameter /3 for k — 30. The upper inset shows the shape 
of the PN barrier before the minimum at f3 = 7.8 (curve 1) 
and after the minimum at /3 = 8.8 (curve 2). The lower inset 
shows more detailed behavior of 5E around its minimum. 

FIG. 17. Normalized velocity of the non-oscillating kink 
motion against the coupling parameter k for (3 = 10. 

FIG. 18. Normalized velocity of the non-oscillating kink 
motion against parameter /3 for k, — 120. 

FIG. 19. Examples of moving antikinks at j3 — 10 and 
k = 100 with velocities: (a) s = 6.227 and (b) s = 2.463. 
Circles show positions of lattice sites. 

FIG. 20. Examples of moving antikinks at k = 120: (a) 
P = 12 and s = 4.154, (b) (i = 15 and s = 4.551, (c) (3 = 20 
and s = 4.734. Circles show positions of lattice sites. 
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This figure "figl2.gif" is available in "gif" format from: 
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